MICROCOPY  RE SOLU l ION  MSI  CHARI 

NA I ION  At  RUM  All  Of  MANUARDS  I'JfcrS  A 


DDC  FILE  COPY  AD  AO  59  650 


r 


AN  ANALYTICAL  STUDY  OF  HEAT  TRANSFER  IN  A TWO-PHASE 
SOLID-LIQUID  MEDIUM  WITH  TIME  VARYING  BOUNDARY  CONDITIONS 
AT  BOTH  THE  MOVING  AND  STATIONARY  BOUNDARIES 


Mark  W.  Nansteel 
Carl  H.  Woigemuth 


V 


Technical  Memorandum 
File  No.  TM  78-175 
June  7,  1978 

Contract  No.  N00017-73-C-1418 

/ 

Copy  No.  (O 


The  Pennsylvania  State  University 
Institute  for  Science  and  Engineering 
APPLIED  RESEARCH  LABORATORY 
P.  0.  BOX  30 

State  College,  PA  16801 


Approved  for  public  release  distribution  unlimited. 


?8 

NAVY  DEPARTMENT 


v*o 


NAVAL  SEA  SYSTEMS  COMMAND 


*9 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whan  Data  F.ntarad) 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


2.  GOVT  ACCESSION  NO.I  3-  RECIPIENT'S  CATALOG  NUMBER 


TM-78-175 


_£N  ANALYTICAL  ^TUDY  OF  J1EAT  .TRANSFER  IN  A JWO- 
_PHASE,  SOLID-LIQUID  J1EDIUM  WITH JTIME JARYING 
BOUNDARY  .CONDITIONS  ’ AT  gOTH  THe’mOVING  AND 


w ■ 0 *1  v'l  i7;i * 3*1 

l|h  I H J I I j I— ■ 


rrnrrl 


Mark  W./Nansteel  j 
Carl  H. /Wolgemuth  / 

9 PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Applied  Research  Laboratory 
P.  0.  Box  30 

State  College,  PA  16801 

It.  CONTROLLING  OFFICE  NAME  ANO  ADDRESS 

Naval  Sea  Systems  Command 
Washington,  DC  20360 


[/ $~)\  NQ&017-73-C-1418 


10.  PROGRAM  ELEMENT,  PROJECT,  TASK 
AREA  & WORK  UNIT  NUMBERS 


(/f),  7 JuHOjfe  WZ  j 


14.  MONITORING  AGENCY  NAME  ft  ADORESS (It  different  from  Controlling  Office)  I 15.  SECURITY  CLASS,  (of  this  report) 


[/£)  J % 


Unclassified 

15a.  (DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


[ 16.  DISTRIBUTION  STATEMENT  (of  IN]  Report) 


Approved  for  public  release;  distribution  unlimited. 

Per  NAVSEA  Sept.  6,  1978. 


17.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  In  Block  20,  If  different  from  Report) 


10.  supplementary  NOTES 


I 19.  KEY  WOROS  (Continue  on  reverse  side  If  necessary  and  Identify  by  block  number) 


heat  transfer 

phase  change 

so  I i d- 1 iqu i d system 


2(JJf  ABSTRACT  (Continue  on  reverse  side  It  necessary  and  Identity  by  block  number) 

The  following  report  provides  a brief  outline  and  description  of  the  steps 
taken  in  formulating  a simple  mathematical  model  for  the  analysis  of  heat 
transfer  in  a pure  substance  undergoing  a change  in  phase  (solid  to  liquid  or 
vice-versa).  The  problem  is  approached  from  the  differential  point  of  view 
and  eventually  the  governing  equations  are  solved  by  a finite  difference 
technique.  The  results  show  fairly  conclusively  that  the  present  method  of 
analysis  is  capable  of  predicting  the  behavior  of  such  a system  to  a 
reasonably  high  degree  of  accuracy. 


DD  1)2""?!  1473  COITION  OF  * NOV  65  IS  OBSOLETE  1 

39100? 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whan  Data  Enlar, 


Subject:  An  Analytical  Study  of  Heat  Transfer  in  a Two-Phase, 

Solid-Liquid  Medium  With  Time  Varying  Boundary  Condi- 
tions at  Both  the  Moving  and  Stationary  Boundaries 

References:  Pages  44  and  45. 

ABSTRACT 


The  following  report  provides  a brief  outline  and  description  of 
the  steps  taken  in  formulating  a simple  mathematical  model  for  the 
analysis  of  heat  transfer  in  a pure  substance  undergoing  a change  in 
phase  (solid  to  liquid  or  vice-versa).  The  problem  is  approached  from 
the  differential  point  of  view  and  eventually  the  governing  equations 
are  solved  by  a finite  difference  technique.  The  results  show  fairly 
conclusively  that  the  present  method  of  analysis  is  capable  of  predicting 
the  behavior  of  such  a system  to  a reasonably  high  degree  of  accuracy. 


^0 


-2- 

TABLE  OF  CONTENTS 

Abstract 

Table  of  Contents  

List  of  Figures  

Nomenclature 

I.  Introduction  

II.  Statement  of  the  Problem  

III.  Basic  Assumptions  in  Model  Formulation.  . . 

IV.  Mathematical  Formulation  

V.  Methods  of  Analysis  

VI.  Nondimensionalization  and  Asymptotic  Solution 

VII.  Finite  Difference  Solution 

VIII.  Results  and  Conclusions 

Appendix 

References 


June  7,  1978 
MWNrlcl 

Page 

....  2 
....  3 

....  4 

....  6 

....  10 
....  13 
....  15 
....  19 
....  22 

....  26 
....  43 
....  44 


-3- 


June  7,  1978 
HWN.-lcl 


LIST  OF  FIGURES 


Figure  Pag 

1.  Sketch  Illustrating  the  Basic  Arrangement 

of  Solid  Layer,  Rigid  Wall,  and  Liquid  Freestream 8 

2.  Variable  Nodal  System 23 

3.  Initial  Condition  for  the  Stefan  Solution, 

No  Solid  Layer. 

Initial  Condition  Specified  in  Part  IV, 

S=S  at  t=0 29 

o 

4.  Illustration  of  Method  for  Comparing  the 

Exact  and  Numerical  Solutions 30 

5.  Comparison  of  Exact  (Stefan)  and  Numerical  Solutions. ...  32 

6.  Asymptotic  Behavior  for  St  *■  0 34 

7.  Comparison  of  Temperature  Gradients  at  the 
Interface,  Actual  Temperature  Profile  and 

Linear  Approximation  35 

8.  Constant  Heat  Flux,  Wall  Temperature  Stepped  Up 37 

9.  Constant  Wall  Temperature,  Heat  Flux  Stepped  Up 38 

10.  Constant  Heat  Flux,  Wall  Temperature  Stepped  Down 39 

11.  Constant  Wall  Temperature,  Heat  Flux  Stepped  Down 40 

12.  Constant  Wall  Temperature,  Step-Wise  Heat 

Flux  Variation 41 


-4- 


June  7,  1978 
MWN : lei 


NOMENCLATURE 


Dimensional  Quantities: 


Symbol 

c 

d 

h 

k 

l 

L 

«c 

S 

S 


Description 
heat  capacity 


Btu 


lbm-  F 


characteristic  length  - ft 

Btu 


convective  conductance  - 


hr-f t -°F 
Btu 


thermal  conductivity  - 0f. 

heat  of  fusion 


Btu 

lbm 


T 

a 

t 


characteristic  length  - ft 

, , ,,  Btu 

convective  heat  flux  - — — 7~7 

hr-f  t 

solid  layer  thickness  - ft 

solid  layer  thickness  at  zero  time  - ft 

temperature  - °F 

melting  temperature  - °F 

freestream  temperature  - °F 

time  - hr 


x 


position  coordinate  measured  parallel  to  the 
solid  wall  - ft 


y 

a 

P 


position  coordinate  measured  perpendicular  to 
the  solid  wall  - ft 

ft2 

thermal  diffusivity  - 

, , lbm 

mass  density  - yjTT 


T 


time  coordinate  for  Stefan  solution  - hr 
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Nondimensional  Quantities: 


Symbol 

b 

Fo 

n 

N 


S* 

St 

Sti 

t* 

6 

n 

e 


Description 

numerical  constant 

Fourier  modulus 

nodal  point  index 

denotes  interface  node 

dimensionless  heat  flux 

dimensionless  heat  flux  for  t<0 

dimensionless  solid  thickness 

Stefan  number 

Stefan  number  for  t<0 

dimensionless  tine 

symbol  denoting  small  magnitude 

similarity  parameter 

dimensionless  temperature 
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I.  INTRODUCTION 

The  study  presented  in  this  paper  was  initiated  in  order  to  gain 
a better  understanding  of  the  phenomina  occurring  in  a two-phase  (solid- 
liquid)  system  under  transient  temperature  conditions.  The  general  phase 
change  problem  considered  here  is  one  of  great  practical  interest  in  Lhe 
areas  of  nuclear  reactor  safety,  the  solidification  of  metal  castings, 
cryogenics,  etc..  Due  to  the  mutual  presence  of  two  phases,  this  problem 
defies  solution  by  techniques  commonly  employed  in  problems  of  conduction 
heat  transfer  because  the  region  in  which  the  solution  is  sought  is 
continuously  changing  with  time. 

The  immediate  application  of  the  present  analysis  concerns  the 
prediction  of  the  behavior  of  the  solid  product  layer  adjacent  to  the 
tube  walls  (heated  side)  in  a Rankine  cycle  steam  boiler-liquid  metal 
reactor.  The  amount  and  rate  of  deposition  of  solid  product  is  sought 
as  a function  of  time.  The  fundamental  goal  or  objective  in  this  study 
then,  is  to  formulate  a simplified  mathematical  model  which  may  be  used 
to  predict  the  rate  of  solidification  or  melting  in  a two-phase  system 
such  as  the  one  described  above  for  a given  set  of  initial  and  boundary 
conditions. 

The  basic  problem  of  heat  transfer  with  phase  change  has  avoided 
general  solution  for  a long  time.  The  analytical  difficulties  encountered 
can  be  mainly  attributed  to  the  nonlinear  character  of  Lhe  governing 
differential  equations.  Numerous  particu 1 a r so  1 ut ions  have  been  presented 
in  the  literature  (see  for  example.  Cars  law  and  Jaeger  (1),  J.  Stefan  (2), 
and  I.  G.  Portnov  (3)  ).  However,  in  general,  the  application  of  these 
solutions  is  limited  to  systems  with  special  boundary  conditions  of  the 
type  which  are  seldom  encountered  in  real-world  melting  and  solidification 
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phenomena . Numerous  approximate  methods  have  also  been  suggested  in- 
cludin'.', variational  techniques  and  integral  methods  similar  to  those 
developed  by  Th.  von  Karraan  for  the  calculation  of  two-dimensional 
laminar  and  turbulent  boundary  layers  in  fluid  dynamics.  The  accuracy 
of  these  approximate  methods  is  in  general  a function  of  the  immediate 
problem  but  they  are  nevertheless  very  powerful.  Relatively  accurate 
results  can  be  secured  with  these  methods  with  the  additional  benefit 
that  the  mathematical  complexity  is  significantly  reduced  in  comparison 
to  the  more  general  solution  techniques.  The  various  methods  of  solu- 
tion will  be  elaborated  upon  later  on  in  this  paper. 

The  brief  discussion  of  solution  methods  given  above  already  implies 
that  in  order  to  obtain  results  for  the  general  melting  and  solidification 
problem  numerical  techniques  must  be  used.  Therefore,  the  major  effort 
in  this  study  will  be  toward  the  development  of  a finite  difference  approx- 
imation for  the  governing  differential  equations.  The  resulting  difference 
equations  will  be  solved  numerically  by  digital  computer.  The  results 
will  be  presented  in  dimensionless  form  so  that  the  effects  of  the 
various  parameters  can  be  more  easily  interpreted.  An  integral  type 
analysis  is  being  developed  at  this  time  which  posseses  some  favorable 
characteristics  even  when  compared  to  the  full  finite  difference  method, 
however,  the  details  of  the  integral  method  will  not  be  discussed 
here. 

II . Statement  of  the  Problem 

Figure  1 shows  a two-phase  system  of  solid  and  liquid  in  which  the 
solid  layer  is  bounded  on  one  side  by  a rigid  wall  and  on  the  other  side 
by  the  liquid  phase.  The  symbols  in  Figure  1 have  the  following  meanings. 
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(LIQUID)  qc  h(T°°  ’ Tm*  INTERFACE 


Figure  1 - Sketch  Illustrating  the  Basic  Arrangement  of 

Solid  Layer,  Rigid  Wall,  and  Liquid  Freestream. 
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h - local  convective  conductance  based  upon  T , the  liquid  freestream 
temperature 

q^  - local  convective  heat  flux  at  the  solid-liquid  interface 
S - solid  layer  thickness 

T - temperature  in  the  solid  layer 

T - equilibrium  fusion  temperature  of  the  substance 

T - local  wall  temperature 

w 

t - time  coordinate  measured  from  some  convenient  origin 
x - position  coordinate  measured  parallel  to  the  rigid  wall 

y - position  coordinate  measured  perpendicular  to  the  wall 

The  wall  is  assumed  to  be  either  a plane  wall  or  one  in  which  the 
local  radius  of  curvature  is  large  with  respect  to  the  local  solid  layer 
thickness,  S.  This  assumption  is  made  because  in  the  present  study  a 
Cartesian  coordinate  system  is  used,  but  the  same  basic  type  of  analysis 
and  procedure  cojld  be  applied  in  either  cylindrical  or  spherical  coordi- 
nate systems.  At  the  solid-liquid  interface,  (y=S),  the  solid  surface 
is  subject  to  the  convective  heat  flux  from  the  warmer  liquid  which  is 
of  course  a function  of  the  local  convective  conductance,  h,  and  the 
freestream  temperature,  T . At  the  wall,  (y=0),  the  layer  is  in  direct 

OO 

contact  with  the  solid  wall  which  is  assumed  always  to  be  at  some  temper- 
ature less  than  the  fusion  temperature,  hence  some  solid  is  always  present. 
The  assumption  of  direct  contact  between  solid  and  wall  implies  the  ab- 
sence of  any  contact  resistance  and  therefore  the  temperature  of  the 
solid  at  y=0  is  equal  to  the  wall  temperature.  It  is  important  to  note 
that  in  this  analysis  no  assumptions  will  be  made  regarding  the  form  of 

the  time  variation  of  the  wall  temperature,  T , and  thr  convective 

w 

heat  flux,  q".  On  the  other  hand  both  T and  q"  must  be  known  as 
V w ^c 
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functions  of  time  before  the  problem  solution  can  begin.  Specification 
of  the  heat  flux  at  the  interface  implies  that  the  fluid  dynamical  problem 
for  the  region  y>S  is  independent  of  the  solid  layer  thickness.  This 
can  usually  be  justified  as  a reasonable  assumption  in  external  flow 
situations  when  S is  small,  but  for  internal  flows  the  variation  in 
solid  thickness  may  have  a significant  effect  on  the  fluid  flow  itself. 

In  such  a case  the  equations  governing  the  flow  in  the  region  y>S  are 
coupled  at  the  solid-liquid  interface  to  the  equations  governing  the 
solid  layer  thickness.  A more  complete  discussion  of  the  assumptions 
made  in  this  analysis  is  given  immediately  below. 


Ill . Basic  Assumptions  in  Model  Formulation 

The  following  simplifying  assumptions  will  apply  throughout  the 
remainder  of  this  paper. 

1.  Heat  transfer  within  the  solid  layer  is  primarily  one-dimensional 
and  in  a direction  normal  to  the  solid  wall.  This  assumption  is 
completely  justified  when  the  solid  layer  thickness,  S,  is  rela- 
tively small  with  respect  to  other  characteristic  dimensions  of 
the  problem.  The  validity  of  this  assumption  for  thin  solid 
layers  is  justified  by  the  following  order  of  magnitude  analysis. 


Consider  the  two-dimensional,  transient  heat  conduction  equation. 


1 = 92T  d2T 

a 3 t c*y'  dx.  2 


III-l 


where  * is  the  thermal  diffusivity  of  the  solid.  Introduce  the 


two  characteristic  length  scales  d and  L where  d is  the  scale 


factor  for  the  coordinate  measured  perpendicular  to  the  wall. 


y.  For  d to  be  a characteristic  length  in  the  transverse  direction 


it  must  be  of  the  same  order  of  magnitude  as  the  solid  layer 
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thickness,  S,  Thus, 

d ~ 0(S). 


L is  the  scale  corresponding  to  the  coordinate  measured  parallel 

to  the  wall,  x.  The  two-dimensional  energy  equation  is  nondimen- 

sionalized  with  the  length  scales  d and  L and  the  fusion  temperature, 

T , yielding, 
m 


8T  d.2 

3t'  'W  '1/  Hx2 


III-2 


where. 


The  scaling  factors  d,  L,  and  T were  chosen  so  that  Y,  X,  and  T 

m 

are  all  of  unity  order. 

Y = J ~ 0(1) 

X = J ~ 0(1) 

T - | ~ 0(1) 

m 

If  the  solid  layer  thickness  is  small  in  comparison  to  its  exten- 
sion in  the  x direction  then. 


r ~ 0(6),  where  6<  < 1. 

L 

The  right  hand  side  of  the  energy  equation  III-2  becomes, 
symbolically. 


32T  . A 2 tfy 

W2  V Tx2 


+ 62 


1 

I2 


l + 6 2 . 
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Therefore  approximately. 


8T 

a t' 


a2,r 

TP 


II-3 


and  the  one  dimensional  assumption  is  seen  to  be  justified  for  thin 
solid  layers.  Note,  however,  that  the  above  argument  breaks  down 
near  edges  or  discontinuities  in  the  wall  because  the  assumption 


d 

L 


< 


< 


1 


is  no  longer  justified  in  these  regions. 

2.  The  wall  temperature,  T , and  the  convective  heat  flux  at 

w 

the  interface,  q",  are  not  dependent  upon  the  solid  layer  thickness, 
c 

The  equation  governing  the  temperature  distribution  in  the  solid 
layer  is  therefore  uncoupled  from  the  energy  equation  governing 
the  wall  temperature  and  the  equations  of  fluid  dynamics  governing 
the  liquid  flow. 

3.  The  material  under  consideration  is  a pure  substance.  The 
complexities  which  arise  due  to  phase  change  temperature  ranges 
and  solid  state  phase  transformations  in  multicomponent  systems 

are  not  considered  here.  With  respect  to  alloy  systems,  for  example, 
Friedman  (4)  has  approached  the  phase  change  problem  successfully 
with  a finite  element  technique,  however,  the  present  analysis  will 
apply  only  to  pure  substances,  homogenious  in  the  solid  phase. 

4.  All  solid  properties  are  temperature  independent.  This 
assumption  will  not  severely  affect  the  accuracy  of  the  results 
in  problems  where  temperature  differences  in  the  solid  are  not 
extreme,  i.e.  for  cases  in  which  there  is  a moderate  degree  of 


subcooling. 
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5.  Thermodynamic  equilibrium  exists  at  the  solid-liquid  inter- 
face, in  other  words,  the  interface  temperature  is  assumed  equal 
to  the  fusion  temperature  at  all  times. 


IV.  Mathematical  Formulation 

For  the  assumptions  stated  above  the  temperature  distribution  in 
the  solid  layer  is  governed  by  the  following  equation. 


The  boundary  conditions  are  obtained  directly  from  Figure  1. 
Initially, 

T = T(y)  , 

S = S . 


IV-1 


0: 


At  the  wall, 
y = 0: 


T = T 


IV-la 

IV-lb 

IV-lc 


At  the  solid-liquid  interface. 


y = S:  T = T . IV-ld 

One  other  condition  must  be  specified  at  the  interface.  This  condition 
relates  the  energy  fluxes  at  y = S to  the  rate  of  latent  heat  evolution 
due  to  phase  change.  In  words:  The  rate  of  energy  transport  by  con- 
vection to  the  interface  from  the  flowing  liquid  plus  the  rate  at  which 
energy  is  evolved  due  to  the  latent  heat  of  phase  change  is  equal  to  the 
rate  of  energy  transport  away  from  the  interface  and  into  the  solid 
layer  by  conduction,  or, 


at  y = S: 


IV- 1 e 
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The  symbols  in  equations  IV-1  have  the  following  meanings: 
a - thermal  diffusivity  of  the  solid 
p - mass  density  of  the  solid 
c - heat  capacity  of  the  solid 
k - thermal  conductivity  of  the  solid 
2.  - latent  heat  of  phase  change 

Equations  IV-1  completely  define  the  present  problem  mathematically. 

No  general  solution  exists  at  this  time  for  the  system,  TV-1. 

The  major  difficulty  with  respect  to  analytical  treatment  arises  from 
the  interface  boundary  condition,  TV-le,  as  shown  below. 

In  general,  the  temperature  in  the  solid  layer  is  both  a function 
of  position  and  time,  thus, 

T = T[y,t] 


and 


at  y = S, 


T = T , and 
m 


Upon  rearrangement  this  becomes » 


dy 

dt 


y=S 


f 8T  1 
9_t 
3T 
9y  J 


dS 

dt 


y=S 


Substituting  this  result  into  condition  IV-le  results  in. 


DT 
9 y ! 


y=S 


q^«  - p«. 


o t 

Jt 

9 y ) 


y*s 
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The  interface  boundary  condition,  IV-le,  is  therefore  nonlinear  and  it 
is  this  nonlinearity  which  causes  analytical  difficulty. 


V.  Methods  of  Analysis 

As  mentioned  earlier,  exact  analytical  solutions  are  available  for 
phase  change  problems  only  of  the  most  fundamental  variety.  These  par- 
ticular solutions  apply  only  for  special  forms  of  conditions  IV-la,  b, 
c,  and  e.  Consider,  for  example,  the  classical  Stefan  solution.  This 
closed  form  solution  applies  to  a semi- inf inite  region  which  initially 
consists  of  a single  phase  (liquid)  at  the  fusion  temperature.  At  zero 

time  the  boundary  temperature,  T is  stepped  down  to  some  lower  temper- 

w 

ature  thus  initiating  formation  of  the  solid  phase.  The  boundary 
temperature  remains  constant  thereafter.  In  this  instance  the  liquid 
phase  is  spatially  isothermal  at  the  fusion  temperature  and  condition 
IV-le  simplifies  to. 


k 


9_T 
3 y 


y=S 


With  this  simplification  and  the  assumption  of  constant  wall  temperature 
a similarity  solution  is  possible  with. 


H = — 1 — • (similarity  variable). 

2/aX 


When  this  relation  is  substituted  into  system  IV-1  an  ordinary  dif- 
ferential equation  results  which  can  be  integrated  to  yield  the  tem- 
perature distribution  in  the  solid  and  the  interface  position. 


T 


T 

w 


erf 


T ) 
w 


2*  at 


erf (b) 


V-l 


S = 2b /at  , 


V-2 
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where  the  constant,  b,  is  given  by, 


b2  c (Tm  - Tw) 

be  erf  (b)  = 

IV  TT 


Obviously,  the  practical  application  of  this  solution  is  severely 
limited  by  the  restrictive  boundary  conditions  which  were  imposed 
in  order  to  yield  the  closed  form  solution,  V-l. 

In  contrast  to  the  Stefan  solution,  I.  G.  Portnov  (3)  obtained 
a series  solution  to  a slightly  more  general  problem  by  Laplace  trans- 
form methods.  However,  according  to  the  paper  of  K.  Stephan  (5)  the 
evaluation  of  even  the  first  few  coefficients  in  the  series  is  very 
cumbersome. 

In  order  to  relax  the  restrictions  imposed  by  an  analytical  solution, 
many  investigations  have  employed  approximate  solution  techniques,  the 
most  popular  of  which  are  the  integral  analyses.  A brief  description 
of  the  integral  method  is  given  below,  however,  a much  more  complete 
account  may  be  found  in  Schlicting  (6)  and  in  (Jzisik  (7). 

The  basic  procedure  used  in  most  integral  analyses  in  heat  transfer 
and  fluid  dynamics  is  to  first  approximate  the  unknown  temperature  or 
velocity  profile  by  some  functional  relationship,  the  form  of  which  is 
chosen  in  order  to  satisfy  the  boundary  conditions  for  the  given  problem. 
The  governing  equation  (equation  of  motion  or  energy)  is  then  integrated, 
using  the  approximate  profile,  with  respect  to  one  of  the  spatial 
variables.  The  resulting  ordinary  differential  equation  is  then  in 
terms  of  the  second  independent  variable  and  may  be  solved  by  standard 
methods.  The  overall  effect  of  the  integral  method  then  is  the  same  as 
the  effect  of  the  similarity  transformation  in  that  the  original  partial 
differential  equation  is  reduced  to  an  ordinary  differential  equation. 
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In  contrast  to  the  similarity  solutions  the  integral  solution  can  only 
be  an  approximation  because  the  assumed  profile  satisfies  the  original 
governing  equation  only  at  certain  specific  points  within  the  given 
region,  at  the  boundaries  for  example.  It  follows  then,  that  the  greater 
the  number  of  boundary  conditions  satisfied  by  the  assumed  profile  the 
more  accurate  the  resulting  integral  solution.  In  many  cases  there  is 
a choice  as  to  which  boundary  conditions  are  to  be  used  in  developing 
the  approximate  profile.  In  this  case  the  conditions  should  be  chosen 
for  the  boundaries  which  are  the  most  critical  in  a given  situation  so 
that  the  assumed  profile  will  be  especially  accurate  in  these  regions. 

In  the  phase  change  problem,  for  example,  it  is  imperative  that  the 
interface  heat  balance  condition,  IV-le,  be  used  informing  the  approx- 
imate profile.  Integral  techniques  have  been  successfully  applied  to 
melting  and  solidification  problems  by  T.  R.  Goodman  (8),  P.  A.  Libby 
and  S.  Chen  (9),  J.  M.  Savino  and  R.  Siegel  (10),  and  K.  Stephan  (5). 

All  of  the  studies  mentioned  here  report  results  which  are  usually 
well  within  the  limits  of  accuracy  normally  required  for  engineering 
purposes.  For  instance,  Stephan  reports  integral  results  which  are 
in  error  by  less  than  2.5  percent  when  compared  to  presumably  exact 
numerical  solutions.  Despite  these  promising  results  the  applica- 
bility of  integral  techniques  has  thus  far  not  been  justified  in  more 
general  phase  change  problems  where  the  temperature  profile  within 
the  solid  may  not  be  approximated  well  by  a simple  polynomial.  Such 
situations  may  arise  when  the  wall  temperature  and/or  the  heat  flux 
are  varied  in  an  irregular  manner.  The  integral  technique  is  a 
powerful  method  of  solution  when  applied  to  phase  change  problems  but, 
more  work  must  be  done  in  order  to  establish  the  limitations  involved 


i 
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when  using  this  type  of  approximate  analysis. 

Another  approximate  method  which  has  been  used  successful ly  in  phase 
change  investigations  involves  the  variational  technique  developed  by 
Biot  (11).  This  method  was  applied  to  the  case  of  constant  wall  temp- 
erature solidification  by  C.  Lapadula  and  W.  K.  Mueller  (12).  The  re- 
sults of  their  work  are  comparable  in  .accuracy  to  the  integral  method 
used  by  Libby  and  Chen  with  the  additional  benefit  that  the  solution  is 
available  in  closed  form.  On  the  whole,  however,  variational  techniques 
are  not  nearly  as  popular  as  integral  methods  and  integral  techniques 
therefore  are  responsible  for  the  bulk  of  the  approximate  solutions. 

The  third  major  category  of  solution  methods  are  the  finite 
difference  techniques.  The  principal  advantage  of  numerical  solution 
is,  of  course,  the  inherent  immunity  of  these  methods  to  nonlinearities, 
time  varying  boundary  conditions,  etc..  But  even  finite  difference 
methods  are  adversely  affected  by  the  moving  boundary  characteristic 
of  all  phase  change  phenomina.  In  fixed  nodal  arrangements,  the  solid- 
liquid  interface  must  move  toward,  reach,  and  pass-by  each  stationary 
temperature  node  so  that  when  the  interface  is  located  between  two 
nodes  its  exact  position  is  unknown  and  therefore  it  must  be  located 
approximately  by  interpolation.  A continuous  account  of  the  fusion 
front  travel  becomes  then  somewhat  of  a problem.  The  difference 
equations  and  interpolation  formulae  can  become  very  complex,  especially 
when  solution  in  two  space  dimensions  is  considered,  for  example,  see 
Springer  and  Olson  (13).  A unique  nodal  scheme  for  one-dimensional 
problems  which  eliminates  the  need  for  interpolation  and  yields  a 
continuous,  accurate  account  of  the  fusion  front  motion  was  suggested 
by  Murray  and  Landis  (14).  In  this  nodal  arrangement,  the  nodal  mesh 
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changes  continuously  with  fusion  front  travel  so  that  the  interface 
position  is  always  coincident  with  a temperature  node,  thus  there  is 
no  need  for  interpolation  and  the  difference  formulae  are  simplified 
considerably. 

Numerical  techniques  provide  a means  of  obtaining  an  accurate 
solution  to  the  phase  change  problem  outlined  in  Part  II  for  realistic 
boundary  conditions.  Because  it  is  a primary  aim  of  this  study  to  pre- 
sent a method  of  solution  which  is  of  a general  nature,  numerical  methods 
will  be  applied  here,  namely,  a finite  difference  technique  similar  to 
that  used  by  Murray  and  Landis. 


VI.  Nondimensionalization  and  Asymptotic  Solution 

The  following  dimensionless  groups  are  defined: 


y*  = *j-  , dimensionless  space  coordinate, 

o 

tot 

t*  = — 5 , dimensionless  time  coordinate  (Fourier  Modulus), 
So 


T - T 

0 * w , dimensionless  temperature. 


g 

S*  = — , dimensionless  solid  layer  thickness, 

o 

where  Sq  denotes  the  solid  layer  thickness  at  zero  time,  t = t*  = 0. 
Substituting  these  dimensionless  variables  into  equations  IV-1  results 
in  the  following  nondimensional  system. 


0<y*<S*: 


0y*2 


I 

Sc  dt* 


[9-1]  + 


39 

3t* 


VI-1 


t*  = 0: 


9 - 0(y*), 

S*  - 1. 


Vl-la 
VI-  lb 


* 


Ad 
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y*  = 0: 


9-0. 


VI-lc 


y*  = S*: 


9 = 1, 


. 99  _ dS* 

St  9y*  q dt*  * 


Vl-ld 


Vl-le 


The  new  parameters  resulting  from  the  nondimensionalization  are 
defined  as. 


j*  = — 1 “ — a , dimensionless  heat  flux, 

1 % p oi 


c “ Tw) 


, Stefan  number. 


Both  q*  and  are  assumed  to  be  time  dependent  in  the  general  case. 

The  steady  state  thickness  is  obtained  from  equations  VI-1  as 
follows.  Under  steady  state  conditions  VI-1  reduces  to. 


Using  VI-lc  and  Vl-ld  in  VI-2, 

9 = s*  (linear  profile).  VI-3 

ss 

Substituting  VI-3  into  the  interface  boundary  condition  Vl-le  results  in. 


tss  * 

S*  “ qss 
ss 


0 or, 


where  q*  , S , and  S*  are  the  values  of  q*,  S . and  S*  respectively 
.ss  tss  ss  t 

under  steady  state  conditions.  When  dimensional  quantities  are  reintro- 
duced, equation  VI-4  is  seen  to  be  merely  a statement  of  the  equality 
of  heat  flux  by  conduction  through  the  solid  layer  to  the  rate  at 
which  energy  is  convected  to  the  interface  from  the  liquid, 

<q")  - k(Tm“Tw) ss  . 

SS  Q 

hSS 
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Probably  one  of  the  most  Important  results  of  the  normalization 
of  equations  IV-1  is  the  appearance  of  the  parameter,  S (termed 
"Stefan  number").  The  physical  interpretation  usually  associated 
with  is  obtained  directly  from  the  definition, 

= c(Tm-Tw)  sensible  heat 

t £ latent  heat 

Thus,  the  Stefan  number  is  often  interpreted  as  a measure  of  the  relative 
effects  of  sensible  and  latent  heat  in  the  solid  layer.  It  is  shown 
by  Lock  (15)  that  in  cases  for  which  convection  effects  are  small  at 
the  interface, 

7T  ~ F (Fourier  Modulus) . 

St  ° 

Therefore  small  values  of  S^  imply  large  values  of  Fq,  which  in  turn 
implies  the  approach  toward  steady  state  conditions  within  the  solid 
layer  (recall  that  a large  Fourier  modulus  indicates  that  heat  capa- 
city effects  are  small  with  respect  to  heat  conduction  effects).  To 
further  illustrate  the  physical  significance  of  this  parameter,  con- 
sider the  case  in  which  the  Stefan  number  assumes  a very  small  value, 
indicating  small  heat  capacity  effects  in  the  solid,  i.e. 


S «1,  and  F »1  . 

t o 

If  convection  effects  are  also  small  and  the  wall  temperature  is  con- 
stant, then,  VI-1  reduces  approximately  to. 


a2e 

9y*  : 


= 0 


VI-5 


Thus,  the  temperature  profile  must  be  approximately  linear, 


0 = 2l 

S* 


VI-6 
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Using  this  result  in  condition  Vl-le  yields, 

dS*  St 
d^  " S*  " q 


VI-7 


Equation  VI-7  is  an  expression  for  S*  as  a function  of  t*  in  the  case 
of  constant  wall  temperature  and  very  small  Stefan  number.  Assuming 
further  a constant  value  of  heat  flux  (q*  = constant)  and  interchanging 
independent  and  dependent  variables  results  in, 


t* 


S* 


iid't) 


St-q*$ 


VI-8 


where  p denotes  a dummy  variable  of  integration. 
Integrating  VI-8, 


f 1 - £1 
St 

_ 

s 

c 


VI-9 


Equation  VI-9  is  a closed  form  expression  for  S*  as  an  implicit  function 
of  t*  for  thp  limiting  case  of  -»  0.  This  solution  is  of  practical 
interest  because  in  many  cases  S is  indeed  very  small  and  VI-9  can 
then  be  expected  to  give  totally  satisfactory  results.  Ice  is  one 
substance  for  which  equation  VI-9  often  applies,  for  ice,  l/c  - 288°F. 


VII.  Finite  Difference  Solution 

In  Figure  2 the  N temperature  nodes  are  spaced  evenly  throughout 
the  solid  layer  with  the  1st  and  Nth  nodes  located  at  the  wall  and  at 
the  solid-liquid  interface  respectively.  In  order  to  maintain  this 
configuration  as  the  interface  moves,  each  temperature  node  (excepting 
the  first)  must  have  a non-zero  velocity  with  respect  to  the  wall  at 


4 


i 


y=0.  The  principal  purpose  for  introducing  the  nodal  mesh,  of  course, 
is  to  allow  the  approximation  of  equations  VI-l  by  finite  difference 
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Figure  2 - Variable  Nodal  System. 
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formulae.  With  the  above  moving  nodal  arrangement,  however,  the  time 
derivative,  as  it  appears  in  equation  VI-1,  cannot  be  written  in  dif- 
ference form  for  the  reasons  outlined  below. 

Writing  the  time  derivative  in  the  more  explicit  form. 


59 

5t* 


d9 

dt* 


DO 

Dt* 


y*  = constant 

emphasizes  the  important  fundamental  difference  between  the  partial 
time  derivative  at  a fixed  point  in  space  and  the  total  time  deriva- 
tive of  a moving  point  or  particle.  In  order  to  express  equation  VI-1 

in  difference  form,  it  must  be  written  in  terms  of  the  total  time 
DO 

derivative,  —7  , which  is  easily  approximated  by  the  nodal  system  in  Figure  2. 
Dt* 

Consider  9 to  be  a function  of  y*  and  t*. 


9 = 9[y*,  t*]  , 


d9 


59 

5 1* 


dt* 


dy* 


d9  = D9_  _ _59__  + dy*  59 
dt*  Dt*  5t*  dt*  5y*  ’ 


or, 


59  = D9_  _ dy*  _59_  . 
5t*  Dt*  dt*  5y* 

Substituting  into  equation  VI-1  results  in. 


VII-1 


D9  m 529  _ 1_  £^t  . . dy*  59 

Dt*  = 5 y*2  “ Sc  dt*  1 dt*  5y* 

Equation  VII-2  is  the  one-dimensional,  transient,  heat 
equation  in  a form  applicable  to  points  moving  in  the  solid 
velocity,  . The  velocity  of  the  nth  node  is  related  to 

face  velocity  by  the  followjng  relationship, 

dy*  _ yjn_  dS* 

dt*  S*  dt*  ’ 

n 


VII-2 

conduction 
layer  with 
the  inter- 

VII-3 


L 
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in  which  y*  denotes  the  position  of  the  nth  node.  Substituting  VII-3 
7n 

into  VII-2  and  rewriting  VII-2  for  the  nth  nodaL  point  yields. 

[9  -i]  + . VII-4 

Dt*  Ty*7  St  dt*  n S*  dt*  dy*  n 


Equation  VII-4  is  substituted  for  equation  VI-1  in  the  original  set 
of  non-dimensional  equations  which  are  now  rewritten  in  a form 
applicable  to  the  variable  difference  mesh  in  Figure  2. 


PQn  329  _!_dSt  Yn  dS*  96 

Dt*  ~ 3y*2  " S,.  dt*  1 n s*  dt*  3y* 

n c 


t*  = 0: 


(•9  = 9 (y*)  , 

^S*  - 1. 


VII-5a 

VII-5b 


VII-5c 


j-Q  = 0N  =1,  vrr“5d 

v*  ~ y*  — S* i 

yn  yN 

'•S 

t 

In  the  numerical  solution  of  equations  VII-5  all  space  derivatives 
were  approximated  with  3-point  central  difference  formulas  with  the 
exception  of  the  gradient  in  the  interface  equation  which  was  approxi- 
mated with  a 3-point  backward  formula.  The  3-point  formula  was  chosen 
in  favor  of  the  simpler  2-point  formula  at  the  interface  in  order  to 
minimize  truncation  error  in  this  critical  region.  Both  time  deriva- 
tives of  temperature  and  solid  thickness  were  approximated  by  2-point 
forward  formulas.  The  detailed  difference  equations  are  listed  in 
the  appendix. 

The  explicit  form  of  the  difference  equations  simplifies  the 
numerical  solution  of  system  VII-5  considerably,  however,  in  an 
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explicit  solution  scheme  the  question  of  stability  must  be  considered. 
According  to  Ozisik  (7)  there  is  no  general  criterion  for  determining 
the  stability  limits  for  nonlinear  problems.  In  most  cases  the  stability 
bounds  are  determined  by  numerical  experiment.  For  the  present  study, 
however,  a somewhat  different  approach  was  taken.  First,  it  was  assumed 
that  the  general  stability  criterion  for  linear,  parabolic,  partial  dif- 
ferential equations  holds  approximately  in  the  nonlinear  case  (see  Ames 
(16)  ).  Then  as  a second  check  on  the  stability  of  the  solution  a pre- 
diction-correction scheme  was  used.  The  corrector  was  applied  to  S*  at 
a given  time  step  until  successive  corrections  differed  by  less  than 
some  predetermined  amount,  or,  if  after  four  corrections  the  above  step 
size  control  criterion  still  had  not  been  satisfied  the  step  size  was 
halved  and  that  iteration  repeated.  In  this  way  a check  on  the  stability 
and  on  the  accuracy  of  the  solution  was  achieved  in  one  step.  Both 
accuracy  and  stability  are  assured  by  the  above  procedure  because  in 
most  cases  accuracy  strongly  implies  stability.  A solution  being 
carried  out  near  the  stability  bound  will  exhibit  a stable  but  slightly 
irregular  behavior  in  advance  of  the  point  where  the  numerical  algor- 
ithm becomes  truly  unstable.  The  above  step  size  control  procedure 
then,  effectively  detects  impending  instability  by  the  slightly  erratic 
behavior  of  the  solution  which  immediately  precedes  it.  This  procedure 
has  been  used  to  obtain  numerical  results  for  a wide  variety  of  initial 
and  boundary  conditions  (including  those  in  Part  VIII)  and  to  date  no 
stability  problems  have  been  encountered. 

V f f I . Resi ilts  and  Conclusions 

The  numerical  results  in  this  section  have  been  included  to 
verify  the  validity  and  accuracy  of  the  finite  difference  technique 
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developed  in  Part  VII  and  to  allow  some  basic  observations  to  be  made 
regarding  the  general  behavior  of  the  physical  system  outlined  in  Part  II. 

All  curves  were  obtained  with  the  moving  nodal  point  network  shown 
in  Figure  2 using  a total  of  ten  temperature  nodes  (eight  internal  nodes). 
This  mesh  size  was  found  to  yield  accurate  results  without  an  excessive 
amount  of  computational  effort.  For  simplicity,  each  of  the  curves 
calculated  (except  for  Figure  5)  begins  with  a solid  layer  initially 
at  steady  state,  i.e., 

dS*  * 

for  t*  < 0,  S*  = 1,  = 0,  and  S^  = q^ 

* 

where  S . and  q.  denote  the  values  of  Stefan  number  and  dimensionless 
ti  ni 

heat  flux  respectively  for  times  less  than  zero.  At  t*  = 0,  either 
S or  q*,  or  both  change  from  their  initial  value  to  some  other  value 
which  may  or  may  not  remain  constant  with  time,  depending  on  the 
specified  boundary  conditions  for  the  immediate  problem. 

The  numerical  solution  can  be  compared  to  the  exact  analytical 
solution  of  Stefan  for  the  special  case  of  q*  = 0.  This  is  a trivial 
example  but  it  nevertheless  provides  an  excellent  means  of  checking  on 
the  performance  of  the  finite  difference  calculation.  It  will  be  re- 
called that  for  the  case  of  zero  solid  thickness  initially  and  liquid 
at  the  fusion  temperature  the  temperature  distribution  in  the  solid 
and  the  interface  position  at  any  time  are  given  by  equation  V-l  and 

V-2  respectively.  r \ 

erf  — X — 


T 


(T 


T ) 
w 


l 2 /co~  ) 

erf (b) 


V-l 


S - 2b  /a  t , 

where  the  constant  b,  is  determined  from, 

.2  b2  ....  S. 

b e erf(b)  = 

/ IT 


V-2 


VIII  -1 
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The  initial  condition  for  the  Stefan  problem  is  unfortunately  inconsis- 
tant  with  the  requirement  of  the  finite  difference  model  that  the  solid 
layer  always  have  a finite  thickness.  The  difference  in  initial  condi- 
tions is  illustrated  in  Figures  3a,  and  3b. 

This  discrepancy  was  accounted  for  in  the  numerical  comparison  by 
allowing  the  Stefan  solution  to  proceed  until  the  solid  layer  thickness, 
S,  calculated  by  Equation  V-2  was  Sq  in  magnitude  at  which  time  the 
numerical  solution  was  begun  using  the  exact  temperature  profile  given 
by  Equation  V-l  as  an  initial  temperature  distribution.  The  procedure 
is  shown  schematically  in  Figure  4. 

In  Figure  4,  T is  the  time  coordinate  corresponding  to  the  Stefan 
solution,  T0  is  the  time  required  for  the  solid  to  build-up  to  a thick- 
ness Sq  as  given  by  Equation  V-2,  and  t is  the  time  coordinate  for  the 
finite  difference  calculation,  the  origin  of  which  coincides  with  x * tc 
From  Figure  4, 


T = t + T. 


so  that  equation  V-2  becomes 


S = 2b/a(t+r0)  . 

However,  by  definition, 


Sq  =■  2b  /ai0  » 


thus. 


To 


4b  2ct 


Substituting  this  into  Equation  VIII-2  , 


4b2  Ci(t  + — — ) , 

4b  2a 


(f  )2  - 4b2  + l . 

O v> 

o o 


VIII-2 
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Figure  3 


Initial  Condition  for 
the  Stefan  Solution, 
No  Solid  Layer. 


Initial  Condition 

Specified  in  Part  IV, 

S=S  at  t=0. 
o 
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(time)  t = 0,  ORIGIN  OF  TIME  FOR 


FINITE  DIFFERENCE 
CALCULATION 


Figure  4 - Illustration  of  Method  for  Comparing  the 
Exact  and  Numerical  Solutions. 


-31- 


June  7,  1978 
MWN : lcl 


and, 


S*2  = 4 b2t*  + 1 . Vlir-j 

Equation  VIII-3  is  the  Stefan  solution  in  terms  of  S*  and  t*,  the 
pertinent  variables  in  the  numerical  solution.  Similarly,  Equation  V-l 
becomes , 

erf 

9 = erf(b)  " 

Equation  VIH-4  was  used  to  generate  the  initial  temperature  profile 
for  the  finite  difference  solution.  For  convenience  the  constant  b was 
set  equal  to  1/2  so  that  Equation  VIt!^3  reduces  to 

S*  = /t*  + 1 , VIII-5 

with  S^  = 0.5923  from  Equation  VIII-1. 

The  results  of  this  comparison  are  shown  in  Figure  5 in  which  the 
data  points  were  plotted  from  the  finite  difference  solution  and  the 
continuous  curve  is  plotted  from  Equation  VIII-5. 

The  discrepency  between  the  two  solutions  is  so  small  that  it 
cannot  be  easily  seen  from  Figure  5.  In  actual lity,  the  error  in  the 
finite  difference  results,  in  the  range  of  time  shown,  is  less  than 
1/10  of  one  percent  based  upon  the  "exact"  results  of  Equation  VIII-5. 
The  same  curve  was  extended  to  a t*  value  of  10,000  to  ascertain  the 
effect  of  large  values  of  S*  (large  increments  between  temperature 
nodes)  on  the  accuracy  of  the  solution.  At  t*  = 10,000,  at  which  time 
the  solid  thickness  is  100  time  its  original  thickness,  the  relative 
error  in  the  numerical  solution  is  still  only  about  l/20th  of  one  per- 
cent. From  the  results  of  this  comparison  it  would  seem  safe  to  assume 
that  the  solution  technique  outlined  in  this  paper  does  provide  results 


2 /t*  + l/4b' 


which  are  highly  reliable. 
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Figure  5 - Comparison  of  Exact  (Stefan)  and  Numerical  Solutions 
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Figure  6 is  a series  of  plots  parametric  in  Stefan  number  verifying 
the  effect  of  very  small  values  of  on  the  solution  for  S*  (see  part 
VI).  For  small  values  of  the  temperature  profile  in  the  solid  layer 
is  virtually  linear  and  the  interface  position  is  predicted  quite  well 
by  Equation  VI-9  denoted  by  the  data  points  in  Figure  6.  However,  as 
the  Stefan  number  is  raised,  the  heat  capacity  in  the  solid  becomes 
more  and  more  important  and  at  S^  values  approaching  unity  the  predic- 
tions of  equation  VI-9,  derived  under  the  assumption  of  a linear  profile, 
become  poor  approximations  at  best. 

Also  note  that  in  each  curve  in  Figure  6,  Equation  VI-9  over- 
predicts the  magnitude  of  the  interface  velocity.  The  reason  for  this 
is  easily  shown  by  the  dimensional  form  of  the  energy  equation. 


92T  _ _1  JT 
2y2  a 3 t 


iv-l 


For  growing  solid  layers,  such  as  the  ones  in  Figure  6,  the  right  hand 
side  of  Equation  IV-1  is  always  negative,  indicating  that  the  temperature 
profile  is  everywhere  concave  down,  see  Figure  7.  The  slope  of  the 
linear  profile  is  seen  to  be  greater  than  that  of  the  actual  profile 
at  the  interface,  and  therefore  substitution  of  the  linear  relation 
VI-6  into  condition  IV-le, 


3 T 
3y 


y=S 


- q"  = 

c dt 


IV-le 


will  result  in  interface  velocities  in  excess  of  those  given  by  the 
finite  difference  solution.  This  qualitative  result  coincides  with 
the  results  of  Figure  6.  Another  very  important  point  concerns  the 
value  of  q*,  the  heat  flux  parameter,  used  in  the  construction  of  the 
curves  in  Figure  6.  If  the  value  of  q*  had  been  chosen  close  to  unity 
very  poor  correlation  would  have  resulted  regardless  of  the  value  of 
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Figure  7 - Comparison  of  Temperature  Gradients  at 
the  Interface,  Actual  Temperature 
Profile  and  Linear  Approximation. 
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Stefan  number  used.  The  relationship  between  and  Fourier  modulus 
cited  in  part  VI  has  no  meaning  unless  convection  effects  at  the  inter- 
face are  sufficiently  small.  This  important  detail  is  implicit  in  the 
development  of  Equation  VI-5,  (see  Reference  (15))  and  has  been  over- 
looked by  at  least  one  author  (Reference  (5)  ). 

Figures  8 through  11  were  included  in  this  section  in  order  to 
give  some  indication  of  how  the  solid  layer  would  behave  when  subjected 
to  simple  step  changes  in  wall  temperature  or  heat  flux.  One  obvious 
conclusion  to  be  made  here  is  that  for  a melting  layer.  Figures  8 
and  9,  an  increase  in  heat  flux  at  the  interface  is  much  more  effective 
in  causing  a response  than  an  increase  in  wall  temperature  (for  the 
same  final  value  of  S*).  This  effect  becomes  more  and  more  pronounced 
as  heat  flux  is  increased  and  S is  decreased. 

In  figure  11,  the  solid  layer  is  growing  (solidification)  and  the 
convective  heat  flux,  q*,  is  seen  to  have  a substantial  effect  on  the 
time  required  to  reach  steady  state  conditions.  As  the  heat  flux  is 
diminished  the  curves  approach  the  characteristic  "square  root"  shape 
as  in  Figure  5.  The  general  tendency  in  all  the  curves  is  toward  short 
response  times  when  the  heat  flux  is  large  and  the  Stefan  number  is  small. 
Figure  12  was  included  to  show  the  variation  in  solid  thickness  for  an 
arbitrary  step-wise  variation  in  convective  heat  flux  which  might  occur 
during  a load  change  in  a device  like  the  boiler-reactor  referred  to 
in  Part  I.  In  this  particular  case  the  wall  temperature  remained 
constant  at  the  initial  value  but  it  could  have  just  as  easily  been  a 
step  or  continuously  varying  function  of  time.  Figure  12  is  there- 
fore an  illustration  of  the  generallity  of  the  solution  technique 


developed  in  Part  VII. 
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Figure  8 - Constant  Heat  Flux,  Wall  Temperature  Stepped  Up. 
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Figure  9 - Constant  Wall  Temperature,  Heat  Flux  Stepped  Up. 
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Figure  12  - Constant  Wall  Temperature,  Step-Wise  Heat  Flux 
Variation. 
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In  conclusion,  the  results  presented  here  show  relatively  con- 
clusively that  the  analysis  developed  in  this  paper  for  dealing  with 
phase  change-moving  boundary  type  problems  is  capable  of  providing  re- 
liable, accurate  data.  This  statement  only  holds  true  if  the  simplifying 
assumptions  set  forth  in  Part  III  can  be  applied  to  the  physical  problem 
without  an  excessive  loss  in  realism.  Probably  the  most  conclusive  way 
to  show  that  the  simplifying  assumptions  do  not  oversimplify  the  actual 
problem  is  by  direct  experiment.  Such  an  experiment  is  at  present  in 
the  developing  stages  and  should  eventually  provide  data  to  directly 
verify  the  analysis  carried  out  here. 
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Difference  Approximations 
Equation  VII-5: 


APPENDIX 


9 r+1  -3r 
n n 

At* 


o: 


3r,  i - 2 3r) 
n+1  n 


dS 


Ay 


k 2 


— [9‘ 

c dt*  1°. 


- 9 


X f 

n yn  dS*(3n+l 
J S*r  dt*  2Ay* 


Equation  VII-5e: 


(3N-2  " + 


2A  y* 


dr,*  m s*r+1  - s*r  . 

dt*  A t* 


In  the  above  difference  formulae  "r"  denotes  the  rth  time  step, 
the  nth  nodal  point,  and  "N"  the  interface  node. 
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